subcortical-cortical interaction FDG fPET
Marina
2025-02-04
#read global mean values from file
library(openxlsx)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.3library(tidyr)
setwd("~/Documents/Promotion/RS-fPET/Analysen/flexible_factorial_multiframes")
#import globals
globals<-c(3319.94816708075,
3418.38236634798,
3497.18808493240,
3577.02786408874,
3659.49101254483,
3756.67744453291,
3809.87836042326,
3898.57802230989,
3999.31684429413,
4067.43533060166,
4141.67033527087,
4227.20607622073,
4302.37920166778,
4387.44581454623,
4470.22750970148,
4523.49026427971,
4592.58739233506,
4638.80380850590,
4657.24528367891,
4721.00248055408,
4752.43926965674,
4810.47495914204,
4811.32067331599,
4859.34592704655,
4898.67618278456,
4943.01107473273,
4949.49331189629,
5007.50625781463,
5060.08885486849,
5074.08189170276,
3463.21449367145,
3555.81843063848,
3639.67267486807,
3734.10450296421,
3801.91657696198,
3895.99314360205,
3988.07394164183,
4041.65026866628,
4110.65100663420,
4192.70609698503,
4292.16732758294,
4377.10338527132,
4501.34535697103,
4567.50772308234,
4678.20219486598,
4740.08440999034,
4811.25442556765,
4950.06230568339,
5021.62014786680,
5072.87333014877,
5129.44649611283,
5167.32522983417,
5239.80294998409,
5266.96495412684,
5340.24395329966,
5332.45039552743,
5368.89692345910,
5442.17776526553,
5495.03385877380,
5531.23103340016,
2724.70573079751,
2801.11826425685,
2841.94192562966,
2895.85016620547,
2979.67118069875,
3037.26791561374,
3084.16229164628,
3141.66176872056,
3213.61968135634,
3251.18932015025,
3354.87844066079,
3415.80844849723,
3442.50561038114,
3522.19001883295,
3616.54687438812,
3650.15476273535,
3676.84572754950,
3726.55949529830,
3758.13322674733,
3766.31175750036,
3811.14747796121,
3847.24413594188,
3875.38767129541,
3916.72257360192,
3935.29158905346,
4008.14502065617,
4007.56326332832,
4017.61149758045,
4064.51578663252,
4080.07384297825,
2757.73491821340,
2843.03943303919,
2901.70930179077,
2947.41837937330,
3023.03942607961,
3082.33762490042,
3153.83327096665,
3196.95147980858,
3247.36186139382,
3296.83185193454,
3335.99736248624,
3439.28392699495,
3502.18362923834,
3559.72237660085,
3601.11486783002,
3671.49462353559,
3753.15602501133,
3816.38627929906,
3895.21495904310,
3943.57474063943,
3989.94693071748,
3997.53670790789,
4018.32091037161,
4053.80621386861,
4087.91674402029,
4123.19040303367,
4147.52051605164,
4187.17629939089,
4203.47891858052,
4236.66892088889,
2764.84129842810,
2831.63049094681,
2868.23098580737,
2943.11071953718,
3006.72909143078,
3050.72280082254,
3128.98465889064,
3191.23539233331,
3225.09571096278,
3298.38408161486,
3361.11676816427,
3431.46253842500,
3487.51032283339,
3553.70501444626,
3626.84549881193,
3684.91259194636,
3729.49641294629,
3819.25979116901,
3857.36478212511,
3901.48035701993,
3944.02609392190,
3979.57068712671,
4011.23892713366,
4037.22939138019,
4072.33272491402,
4098.10871243931,
4119.79011274126,
4136.17873738423,
4179.31292774336,
4202.42435953202,
2723.64144111185,
2773.95236388179,
2833.74471530761,
2898.16178133546,
2949.09968331363,
2992.62720962392,
3063.71490928920,
3109.14662896132,
3196.59198556991,
3250.88092906892,
3313.01084222412,
3375.63362340496,
3447.15665147918,
3499.09516315241,
3569.49313945483,
3629.23166233416,
3692.81618154316,
3751.45473922890,
3838.66954179513,
3866.62543219542,
3933.03838373645,
4022.04457899870,
4092.28405984656,
4147.63029732313,
4216.44309019632,
4260.31808141891,
4302.52608464204,
4308.82993789807,
4370.55669470633,
4422.24245858151,
4579.32950466242,
4695.54596611744,
4826.28313527527,
4941.78658796845,
5070.41525124966,
5204.68016355390,
5312.88557906887,
5414.32055575324,
5552.69664562181,
5665.88084524768,
5770.85689811420,
5910.12266769570,
6055.05070816002,
6183.28227574275,
6321.77272297455,
6453.90573384614,
6567.16347121194,
6693.51714455732,
6845.42074371483,
6978.83467842561,
7147.19165375549,
7267.04924369983,
7397.61052857497,
7550.31721321434,
7682.23140108729,
7817.06622310416,
7914.70280347762,
7989.47077624116,
8079.99002331678,
8132.24249384481,
3056.76942796260,
3079.33489472465,
3069.81469325774,
3102.69448316216,
3114.07162722644,
3123.54793197637,
3130.63511890437,
3131.36310117853,
3162.91902240081,
3145.20460355814,
3165.47996677631,
3180.30155220493,
3193.02739021486,
3193.88239047159,
3214.11538545042,
3210.46038037488,
3232.77169162684,
3235.82289734160,
3252.38705476565,
3260.45284295437,
3259.55073254645,
3272.28148583680,
3290.35088447513,
3301.16464624867,
3309.35639926200,
3309.53466245413,
3318.78537516829,
3334.99442146378,
3334.28818751421,
3350.24857297137,
2700.62312867496,
2755.06339996951,
2834.28807717885,
2884.48202811316,
2979.37882914823,
3031.58182553708,
3102.96287881290,
3159.86628633659,
3212.49067008733,
3313.03134998421,
3393.99219187691,
3456.27230331572,
3533.04683923934,
3603.42854410941,
3673.28194057120,
3767.63311059231,
3850.06582376081,
3913.29520744008,
4010.54489933994,
4080.26492847487,
4159.40926760452,
4242.93840701548,
4304.85606515297,
4396.46639206363,
4490.70490701455,
4563.65927237689,
4626.59600565548,
4727.00286163860,
4793.77030095492,
4881.22658436441,
4672.95514577808,
4783.40033970401,
4886.54963545112,
5000.74606887873,
5136.50447848172,
5256.37723139465,
5356.75318282963,
5486.23087647404,
5597.94950718557,
5734.09424436729,
5831.52298771340,
5945.27164279603,
6070.81041059454,
6180.51836916879,
6320.35049299492,
6444.03620962951,
6532.63823219008,
6674.69319906543,
6806.41546297530,
6870.46406698631,
7032.25620628118,
7172.34339203785,
7273.25446131122,
7401.24180403076,
7502.58459165665,
7684.58555921490,
7775.45281498599,
7920.94821575943,
8049.51885025144,
8185.58483663639,
2095.06239833971,
2170.62447948463,
2197.44574060888,
2224.80754885885,
2284.58267989466,
2340.99848651292,
2390.47434655120,
2439.74748694997,
2454.07915972278,
2492.87999071611,
2484.68442805821,
2582.48439599537,
2644.65058151161,
2713.51234752882,
2752.70857391481,
2835.28138078382,
2849.18176651456,
2820.99475163803,
2885.57492724856,
2980.16734319188,
3040.64176818402,
3096.62513823249,
3131.98450663043,
3186.28623348476,
3227.99272130983,
3282.51509868892,
3398.71277817170,
3425.58471209891,
3432.67547188267,
3551.87896379541,
3257.50458367897,
3317.40485680203,
3365.07081818868,
3447.94980640515,
3549.27168614472,
3596.00705635560,
3717.30031242694,
3802.13508250150,
3917.76123578062,
3960.42220524238,
4063.21324091555,
4131.31270674617,
4162.88910898933,
4202.50848785036,
4252.12673433285,
4287.91308468453,
4312.40377621723,
4368.89807513131,
4384.84724778483,
4432.64163500377,
4504.28195761995,
4528.48860036807,
4587.52904059029,
4665.55910797011,
4770.64986654619,
4846.35881427521,
4930.32911283703,
5039.10203474194,
5049.84798240388,
5087.07554471904,
3100.84925019087,
3179.93087804215,
3254.30661126453,
3336.54711194584,
3412.62460927669,
3495.78105421507,
3579.98975234024,
3654.23387396665,
3726.73680647839,
3817.79956559307,
3900.39001548916,
3978.66007161610,
4058.09336105587,
4155.00022625346,
4229.02555258606,
4327.28173478395,
4398.08719229951,
4497.73414778229,
4572.18681449414,
4660.38519731358,
4769.72161316680,
4852.92833046445,
4945.37497241966,
5019.09043501711,
5105.12595091168,
5191.87801954082,
5294.72453935502,
5418.08700713833,
5438.25379113736,
5506.43308116922,
2479.84633469280,
2558.45863745364,
2586.74015175475,
2640.22563135263,
2677.61654918587,
2684.53495994127,
2761.02425179149,
2829.77632823087,
2894.06121548796,
2959.52543335916,
3002.91526262532,
3065.29641932065,
3140.64704133912,
3185.89986248209,
3255.77755388229,
3297.32204327621,
3398.19865565633,
3436.13304742216,
3498.99201231784,
3579.74691009578,
3675.43046536305,
3747.11089701406,
3818.35303221822,
3894.58254220171,
3872.62796097249,
3916.38682537284,
3958.77267715184,
4041.37263780339,
4102.37920044991,
4136.85644470996,
4441.95439767401,
4534.30210515798,
4678.90025586019,
4767.49654857554,
4890.32460045874,
5003.42653992291,
5134.03690854965,
5256.05442907187,
5376.58332146329,
5497.88288118729,
5618.94434752291,
5732.84776340330,
5848.96387339833,
5980.89973344712,
6083.74199874933,
6211.08955242319,
6346.15326389427,
6474.02361414170,
6578.85093397235,
6703.96094708056,
6861.13603246712,
6953.17563287754,
7093.53304435392,
7205.89217918149,
7344.35016682937,
7471.63393466732,
7590.95893123734,
7729.12903252164,
7857.34999368733,
7943.76762260336,
2845.54900754000,
2917.05861909940,
2987.82080796886,
3088.78939026572,
3166.29768659137,
3294.31516871021,
3364.41859549549,
3424.46823041624,
3518.04184246446,
3589.37473511420,
3671.54886693156,
3708.77400784638,
3854.34305464441,
3955.07087503823,
4015.42241716590,
4093.24139781126,
4186.00609458063,
4247.02843295876,
4356.07447440638,
4474.60110418588,
4524.00969304423,
4610.55937407860,
4690.56816891014,
4769.46563289049,
4883.92801080927,
4941.39255318814,
5045.64882476297,
5139.80823816589,
5209.26695196210,
5279.03926299111,
4292.51739610408,
4389.49590255018,
4501.81899080335,
4618.92418680273,
4723.12900484050,
4814.08011313969,
4928.43685016322,
5009.06197088842,
5109.35900078647,
5218.30248947565,
5321.69968506417,
5433.34009990898,
5555.47059522491,
5659.89446570111,
5753.96834632869,
5872.75166724424,
5963.51080422932,
6067.71432420295,
6174.02272654688,
6282.47734960793,
6391.03490259500,
6478.86162117558,
6629.30818143491,
6722.97071202891,
6826.83871857741,
6939.83085226035,
6982.91069179954,
7068.08218962183,
7151.68976043225,
7193.55096766104,
5887.87741559224,
6056.86700499497,
6197.39648564736,
6347.75534685844,
6469.03344320978,
6626.83431506232,
6746.21022048908,
6911.43442595896,
7068.51452351016,
7213.45139808186,
7386.82150453702,
7515.36539326744,
7662.24517212888,
7849.03950296941,
7979.44192855476,
8178.69989581313,
8344.61598216003,
8472.14063858022,
8647.75666381938,
8777.26311492773,
8961.04365313878,
9099.22656682185,
9280.20790468723,
9435.01200524830,
9615.09029284948,
9754.41708808525,
9918.30715733493,
10087.8871693458,
10233.5587514870,
10396.0939874438,
4579.33252094656,
4695.59838833431,
4826.31285687290,
4941.81619588835,
5070.39000863408,
5204.71286680576,
5312.91851661789,
5414.29327049736,
5552.70054926923,
5665.85228663364,
5770.89415883801,
5910.12671945010,
6054.95325438077,
6183.35557584894,
6321.77599960135,
6453.94710444522,
6567.24369094516,
6693.59671744208,
6845.39040265081,
6978.79944364394,
7147.15838049054,
7267.01878559355,
7397.61876803269,
7550.28485210578,
7682.28374193943,
7817.02954160814,
7914.70926271497,
7989.47730817632,
8079.95198319183,
8132.29674717295,
2616.55968027273,
2698.68085193695,
2744.80825447971,
2807.21240724668,
2869.56236155191,
2935.36956747497,
3003.76290041638,
3060.36600456702,
3146.01010207412,
3201.57355448887,
3253.39179866598,
3340.17771251816,
3401.46047397038,
3478.61845134193,
3546.64775305503,
3607.64234635013,
3642.99560786268,
3729.26997034677,
3808.84973295660,
3841.80296899294,
3942.10504904216,
4000.50626008020,
4054.55306936458,
4128.74750839528,
4196.84634100951,
4245.88514064596,
4305.38652995404,
4385.65969013307,
4417.00809879208,
4478.64807433623,
2910.08967780965,
2990.63038056923,
3042.83788963236,
3128.31449581930,
3199.56234745909,
3279.30696810227,
3365.80922260096,
3410.87173234437,
3513.94324672892,
3574.90224974548,
3636.05443030603,
3721.34988897654,
3789.69704889620,
3853.17000109644,
3946.32729302398,
4008.27140365347,
4093.37674560777,
4171.94404835394,
4234.02814473147,
4322.59180003881,
4421.90279467229,
4486.63703164228,
4562.90335489862,
4672.59998682407,
4731.36916801668,
4822.48990138430,
4883.83057986733,
5004.51795504542,
5071.95059849855,
5165.68930788593,
1595.02197722526,
1638.82506120646,
1662.39544284796,
1704.78425724697,
1740.44876228869,
1780.07226154444,
1813.72924692219,
1856.26259312770,
1906.68252340686,
1932.43421617007,
1967.22414590665,
2030.81182339258,
2053.92154281495,
2101.58751286084,
2136.26950565587,
2185.34959929394,
2239.77938375152,
2238.16026894445,
2290.93066517181,
2336.21477256619,
2370.54632018562,
2431.34782680987,
2461.42970551090,
2463.31816783571,
2526.64268090659,
2576.27909282674,
2606.27789975522,
2657.63074971418,
2710.45023055670,
2728.84888161173,
2830.38874262949,
2922.19570724204,
2972.47775759304,
3052.22312356074,
3091.90490437638,
3169.83574380782,
3230.84049794811,
3281.36689162791,
3352.91744841731,
3407.36221083120,
3481.71626080718,
3534.25284706631,
3589.16389475647,
3664.05356446221,
3714.72174200627,
3770.83051925415,
3859.09081227308,
3950.61670148412,
4025.95039420709,
4099.30461256301,
4198.73201656123,
4264.35058088041,
4334.36027917912,
4416.65824384271,
4456.74174673857,
4532.41071578495,
4636.55247795948,
4683.28337572792,
4751.12525029651,
4786.72295484841,
2616.90226029821,
2687.99153768661,
2765.20163477295,
2825.36907028436,
2888.27903816430,
2952.15180143661,
3027.93972279379,
3097.90551229701,
3180.62563989970,
3241.68459941748,
3310.43813694709,
3387.71914287298,
3429.33094047800,
3527.34381768486,
3600.61469575232,
3658.24601944496,
3740.16446202610,
3805.26611363195,
3867.99829405987,
3919.16869566141,
4012.12919413232,
4075.38318184229,
4153.24682474383,
4209.56464903985,
4281.06435434925,
4368.61650095518,
4441.33064997784,
4511.04393370240,
4591.35236671113,
4660.41253262556,
3401.74954024533,
3470.92629357341,
3546.04232399719,
3658.75923587124,
3728.84040983163,
3808.00988761300,
3908.16113098115,
3988.31563402157,
4086.06295940815,
4178.56025872323,
4263.53224166557,
4342.00505681526,
4433.13288197271,
4512.91342570933,
4600.93605821228,
4671.44885840151,
4797.44882853567,
4873.32939619339,
4984.70755605811,
5051.13892464322,
5122.22761752905,
5222.70992542209,
5300.63075602789,
5428.82151907296,
5528.33327599924,
5595.13578139300,
5689.25021625255,
5766.28779112918,
5861.11393302149,
5926.53774068577,
3222.59646860021,
3338.35298544917,
3408.64082530093,
3475.53994809764,
3596.05833640940,
3671.70707456757,
3744.54674953722,
3857.89743163596,
3961.39389545491,
4050.15292937940,
4132.45286845147,
4201.15240970621,
4276.25891870123,
4352.26611336383,
4477.30659872174,
4564.13647048005,
4648.82179053384,
4707.84952303158,
4827.00996343507,
4888.82586608797,
4989.13281016028,
5052.23910247192,
5172.91218690437,
5269.21922942399,
5346.00625930406,
5436.43195410137,
5545.65061055734,
5595.30060931343,
5700.27646432542,
5765.46523956464,
2457.66540424913,
2526.78924200534,
2597.57297735770,
2651.42326132279,
2709.52522562361,
2764.50236610288,
2828.78536941098,
2887.30069481355,
2955.35541733852,
3042.40184578059,
3108.94890230372,
3170.88721522308,
3239.38520674147,
3312.77400626023,
3385.02818353278,
3486.98273682513,
3567.72053557674,
3655.75021520328,
3748.92138712446,
3753.46501947168,
3854.78289984392,
3875.25529380633,
3971.53510476111,
4031.87521561430,
4103.70865248656,
4170.53428709115,
4233.28388026583,
4309.40304256763,
4398.00708602172,
4406.17706750327)
#read data
FlexData<-read.xlsx(xlsxFile = "flex_fact_new2025_02_05.xlsx", sheet = "Tabelle1",
colNames = TRUE)
#devide values element-wise by global mean
gmFlexData<-FlexData/globals
#define subject variable subject id
subject_id<-c(rep("RS_fPET001",30), rep("RS_fPET002",30), rep("RS_fPET003",30),rep("RS_fPET004",30), rep("RS_fPET005",30),rep("RS_fPET006",30),rep("RS_fPET007",30),rep("RS_fPET008",30),rep("RS_fPET010",30),
rep("RS_fPET011",30),rep("RS_fPET012",30),rep("RS_fPET014",30),rep("RS_fPET015",30),rep("RS_fPET016",30),
rep("RS_fPET017",30),
rep("RS_fPET018",30),rep("RS_fPET019",30),rep("RS_fPET020",30), rep("RS_fPET021",30),rep("RS_fPET022",30),rep("RS_fPET023",30),rep("RS_fPET024",30),rep("RS_fPET025",30),
rep("RS_fPET026",30),rep("RS_fPET027",30),rep("RS_fPET028",30),rep("RS_fPET029",30))
#add and convert record id into factor
gmFlexData$record_id<-subject_id
gmFlexData$record_id<-factor(gmFlexData$record_id)
#add diagnose to_id gm removed data set
gmFlexData$diagnose[gmFlexData$record_id=="RS_fPET001"|gmFlexData$record_id=="RS_fPET002"|gmFlexData$record_id=="RS_fPET003"|gmFlexData$record_id=="RS_fPET005"|gmFlexData$record_id=="RS_fPET007"|gmFlexData$record_id=="RS_fPET008"|gmFlexData$record_id=="RS_fPET010"|gmFlexData$record_id=="RS_fPET014"|gmFlexData$record_id=="RS_fPET016"|gmFlexData$record_id=="RS_fPET022"|gmFlexData$record_id=="RS_fPET025"|gmFlexData$record_id=="RS_fPET026"|gmFlexData$record_id=="RS_fPET028"|gmFlexData$record_id=="RS_fPET029"]= "PD"
gmFlexData$diagnose[gmFlexData$record_id=="RS_fPET004"|gmFlexData$record_id=="RS_fPET006"|gmFlexData$record_id=="RS_fPET011"|gmFlexData$record_id=="RS_fPET012"|gmFlexData$record_id=="RS_fPET015"|gmFlexData$record_id=="RS_fPET017"|gmFlexData$record_id=="RS_fPET018"|gmFlexData$record_id=="RS_fPET019"|gmFlexData$record_id=="RS_fPET020"|gmFlexData$record_id=="RS_fPET021"|gmFlexData$record_id=="RS_fPET023"|gmFlexData$record_id=="RS_fPET024"|gmFlexData$record_id=="RS_fPET027"]= "HC"
#convert diagnose into factor
gmFlexData$diagnose<-factor(gmFlexData$diagnose)
summary(gmFlexData$diagnose)## HC PD
## 390 420#add timepoints
time<-c(rep(61:90, 27)) #repeat time points 61- 90 27 x
gmFlexData$time<-time#save activity-time plots for all subjects with mean line per group
for (i in 1:ncol(gmFlexData[,1:12])) {
MeanTime <- aggregate(gmFlexData[,i] ~ time + diagnose, data = gmFlexData, FUN = mean)
MeanTime_df <- cbind(MeanTime[-ncol(MeanTime)], MeanTime[[ncol(MeanTime)]])
colnames(MeanTime_df)<-c("time","diagnose",colnames(gmFlexData)[i])
for (j in 3:ncol(MeanTime_df)){
plot <- ggplot(gmFlexData, aes(x = time, y = gmFlexData[,i])) +
geom_line(aes(color = diagnose, group = record_id), linetype = 3) +
geom_line(data = MeanTime_df, aes(x = time, y = MeanTime_df[,j], color = diagnose), linewidth = 1) +
theme_bw()+
labs(title = "", x = "Time (Min)",
y = expression("Relative ["^18*"F]-FDG Uptake"))+
scale_color_manual(values=c("#1F78B4", "#BFD200"))+
theme(legend.position = c(0.1,0.85))+
theme(legend.title = element_blank())+
theme(text = element_text(size = 25),
axis.text.x = element_text(size=20),
axis.text.y = element_text(size=20))+
show(plot)
ggsave(plot,file=paste0("gmFlexData",colnames(gmFlexData)[i],".png"))}}## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## function (x, y, ...)
## UseMethod("plot")
## <bytecode: 0x127c88c40>
## <environment: namespace:base>## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
setwd(“~/Documents/Promotion/RS-fPET”)
#change wd
#mean uptake per subject
MeanUptake <- aggregate(gmFlexData[,c(1:12)], by= list(subject=gmFlexData$record_id), FUN = mean)
#add diagnose to_id gm removed data set
MeanUptake$diagnose[MeanUptake$subject=="RS_fPET001"|MeanUptake$subject=="RS_fPET002"|MeanUptake$subject=="RS_fPET003"|MeanUptake$subject=="RS_fPET005"|MeanUptake$subject=="RS_fPET007"|MeanUptake$subject=="RS_fPET008"|MeanUptake$subject=="RS_fPET010"|MeanUptake$subject=="RS_fPET014"|MeanUptake$subject=="RS_fPET016"|MeanUptake$subject=="RS_fPET022"|MeanUptake$subject=="RS_fPET025"|MeanUptake$subject=="RS_fPET026"|MeanUptake$subject=="RS_fPET028"|MeanUptake$subject=="RS_fPET029"]= "PD"
MeanUptake$diagnose[MeanUptake$subject=="RS_fPET004"|MeanUptake$subject=="RS_fPET006"|MeanUptake$subject=="RS_fPET011"|MeanUptake$subject=="RS_fPET012"|MeanUptake$subject=="RS_fPET015"|MeanUptake$subject=="RS_fPET017"|MeanUptake$subject=="RS_fPET018"|MeanUptake$subject=="RS_fPET019"|MeanUptake$subject=="RS_fPET020"|MeanUptake$subject=="RS_fPET021"|MeanUptake$subject=="RS_fPET023"|MeanUptake$subject=="RS_fPET024"|MeanUptake$subject=="RS_fPET027"]= "HC"
for (i in 2:ncol(MeanUptake[,2:14])){
boxplot(MeanUptake[,i]~diagnose, data=MeanUptake, main = paste0("Mean",colnames(MeanUptake)[i]))
results<-t.test(MeanUptake[,i]~diagnose, data=MeanUptake)
print(results)}
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = 3.4257, df = 22.716, p-value = 0.002339
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## 0.04290399 0.17393796
## sample estimates:
## mean in group HC mean in group PD
## 1.492089 1.383668
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = 2.2811, df = 21.302, p-value = 0.03293
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## 0.01016617 0.21800837
## sample estimates:
## mean in group HC mean in group PD
## 2.061555 1.947468
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = 2.8208, df = 22.517, p-value = 0.009816
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## 0.03491988 0.22785672
## sample estimates:
## mean in group HC mean in group PD
## 1.736867 1.605478
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = 3.0073, df = 22.5, p-value = 0.00638
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## 0.04186964 0.22715419
## sample estimates:
## mean in group HC mean in group PD
## 2.071889 1.937378
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = 3.2926, df = 24.783, p-value = 0.002981
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## 0.05223899 0.22694852
## sample estimates:
## mean in group HC mean in group PD
## 2.007700 1.868106
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = 1.6719, df = 18.286, p-value = 0.1116
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## -0.03026672 0.26745263
## sample estimates:
## mean in group HC mean in group PD
## 2.147968 2.029375
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = 2.6173, df = 21.182, p-value = 0.01603
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## 0.02691038 0.23456181
## sample estimates:
## mean in group HC mean in group PD
## 1.861839 1.731103
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = 3.1029, df = 23.36, p-value = 0.004954
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## 0.03958414 0.19752980
## sample estimates:
## mean in group HC mean in group PD
## 2.051667 1.933110
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = 3.2518, df = 23.352, p-value = 0.003468
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## 0.03718395 0.16691831
## sample estimates:
## mean in group HC mean in group PD
## 1.517768 1.415717
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = -2.426, df = 17.77, p-value = 0.02615
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## -0.25449267 -0.01815847
## sample estimates:
## mean in group HC mean in group PD
## 2.264076 2.400402
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = -3.4362, df = 22.665, p-value = 0.002286
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## -0.26638508 -0.06607342
## sample estimates:
## mean in group HC mean in group PD
## 1.909971 2.076201
##
## Welch Two Sample t-test
##
## data: MeanUptake[, i] by diagnose
## t = -2.8365, df = 24.157, p-value = 0.009087
## alternative hypothesis: true difference in means between group HC and group PD is not equal to 0
## 95 percent confidence interval:
## -0.19255678 -0.03038962
## sample estimates:
## mean in group HC mean in group PD
## 2.243670 2.355143#association subcortical vs cortical left SN vs. cortical hypermetbolism precentral right
plot(gmFlexData$HC_gr_PD_m10m16m16,gmFlexData$HC_kl_PD_p5m28p61, col=gmFlexData$diagnose)
#cave: all data points
model<-lm(gmFlexData$HC_kl_PD_p5m28p61 ~ gmFlexData$HC_gr_PD_m10m16m16)
abline(model)
ggsave("association_subcorticalvscortical.png")## Saving 7 x 5 in imagesummary(model)##
## Call:
## lm(formula = gmFlexData$HC_kl_PD_p5m28p61 ~ gmFlexData$HC_gr_PD_m10m16m16)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38537 -0.09142 -0.01741 0.06612 0.45476
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.41846 0.04547 53.182 <2e-16 ***
## gmFlexData$HC_gr_PD_m10m16m16 -0.29410 0.03147 -9.345 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1446 on 808 degrees of freedom
## Multiple R-squared: 0.09754, Adjusted R-squared: 0.09642
## F-statistic: 87.33 on 1 and 808 DF, p-value: < 2.2e-16#corr whole group mean values
cor.test(MeanUptake$HC_kl_PD_p5m28p61,MeanUptake$HC_gr_PD_m10m16m16)##
## Pearson's product-moment correlation
##
## data: MeanUptake$HC_kl_PD_p5m28p61 and MeanUptake$HC_gr_PD_m10m16m16
## t = -3.3098, df = 25, p-value = 0.002836
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.7703958 -0.2176128
## sample estimates:
## cor
## -0.5519803#with ggplot
# Change point shapes and colors
ggplot(data=gmFlexData, aes(x=HC_gr_PD_m10m16m16, y=HC_kl_PD_p5m28p61)) +
geom_point(aes(color=diagnose),size=3,alpha=0.5)+
geom_point(data = MeanUptake,
aes(x = HC_gr_PD_m10m16m16, y=HC_kl_PD_p5m28p61, fill=diagnose),
shape=21, size=3, color="black")+
geom_smooth(data = MeanUptake,
method='lm', formula= y~x, color="black")+
scale_color_manual(values=c("#1F78B4", "#BFD200"))+
scale_fill_manual(values=c("#1F78B4", "#BFD200"))+
theme_bw()+
theme(legend.title = element_blank())+
labs(title = "", x = "left SN",
y = "Motor cortex")+
theme(text = element_text(size = 30),
axis.text.x = element_text(size=20),
axis.text.y = element_text(size=20))+
theme(legend.position = c(0.9, 0.8), legend.title=element_blank())+
geom_text(x=1.6, y=2.47, label="r = -0.55, p = 0.003", color="black", size=10 )
ggsave("association_subcorticalvscortical_SN_M1.png")## Saving 7 x 5 in image#ROC curve
library(Epi)## Warning: package 'Epi' was built under R version 4.3.3MeanUptake$Fall<-MeanUptake$diagnose=="PD"
MeanUptake$Fall<-factor(MeanUptake$Fall)
str(MeanUptake)## 'data.frame': 27 obs. of 15 variables:
## $ subject : Factor w/ 27 levels "RS_fPET001","RS_fPET002",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ HC_gr_PD_m10m16m16 : num 1.38 1.36 1.24 1.46 1.4 ...
## $ HC_gr_PD_m14m98m4 : num 1.95 1.82 1.99 2.09 1.77 ...
## $ HC_gr_PD_m16_m2_p12: num 1.78 1.7 1.28 1.74 1.63 ...
## $ HC_gr_PD_m34_p48_p2: num 1.98 1.97 1.83 1.95 2.12 ...
## $ HC_gr_PD_m52m58p24 : num 1.97 1.84 1.77 1.98 1.69 ...
## $ HC_gr_PD_p14m96m16 : num 2.1 1.83 2 2.26 1.77 ...
## $ HC_gr_PD_p22p2p10 : num 1.86 1.64 1.69 1.81 1.87 ...
## $ HC_gr_PD_p46m62p36 : num 2.01 1.88 2 2.11 1.77 ...
## $ HC_gr_PD_p8m16m16 : num 1.48 1.53 1.23 1.52 1.43 ...
## $ HC_kl_PD_m33m9p4 : num 2.37 2.4 2.34 2.18 2.53 ...
## $ HC_kl_PD_p5m28p61 : num 2.21 1.99 2.4 1.98 1.92 ...
## $ HC_kl_PD_p35m2p3 : num 2.36 2.34 2.45 2.2 2.53 ...
## $ diagnose : chr "PD" "PD" "PD" "HC" ...
## $ Fall : Factor w/ 2 levels "FALSE","TRUE": 2 2 2 1 2 1 2 2 2 1 ...ResultsROC<-with(MeanUptake,
ROC(test=HC_kl_PD_p5m28p61, stat=Fall, plot="ROC",MI=F))
MeanUptake$HC_kl_PD_p5m28p61## [1] 2.209212 1.990884 2.401726 1.979455 1.918860 1.718265 2.010278 2.182069
## [9] 2.128547 1.962481 1.963996 2.270269 1.839517 2.045838 2.003197 1.999362
## [17] 1.921927 1.897130 2.010274 1.994589 1.869792 1.723270 1.874417 2.054764
## [25] 1.940962 2.080473 1.904883MeanUptake$Fall<-as.numeric(as.factor(MeanUptake$Fall))-1
#ROC curve with multiple variables
#all clusters
logit_model<-glm(Fall ~ HC_kl_PD_p5m28p61 + HC_gr_PD_p8m16m16+HC_gr_PD_m34_p48_p2+HC_kl_PD_m33m9p4+HC_gr_PD_m16_m2_p12+ HC_gr_PD_m10m16m16+HC_gr_PD_p14m96m16+HC_gr_PD_p46m62p36+HC_kl_PD_p35m2p3+HC_gr_PD_m14m98m4+HC_gr_PD_m52m58p24+HC_gr_PD_p22p2p10, data=MeanUptake)
MeanUptake$pred_probs<-predict(logit_model,type="response")
ResultsROC<-with(MeanUptake,
ROC(test=pred_probs, stat=Fall, plot="ROC",MI=F))## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
library(pROC)## Warning: package 'pROC' was built under R version 4.3.3## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
##
## The following objects are masked from 'package:stats':
##
## cov, smooth, var#define object to plot
rocobj <- roc(MeanUptake$Fall,MeanUptake$pred_probs)## Setting levels: control = 0, case = 1
## Setting direction: controls < casesauc <- round(auc(MeanUptake$Fall, MeanUptake$pred_probs),4) ## Setting levels: control = 0, case = 1
## Setting direction: controls < cases#create ROC plot
ggroc(rocobj, colour = '#00BFC4', size = 2) +
ggtitle(paste0('ROC Curve ', '(AUC = ', auc, ')'))+
theme_minimal()
#reduced clusters
logit_model_hypo<-glm(Fall ~ HC_gr_PD_p8m16m16+HC_gr_PD_m34_p48_p2+HC_gr_PD_m16_m2_p12+ HC_gr_PD_m10m16m16+HC_gr_PD_p14m96m16+HC_gr_PD_p46m62p36+HC_gr_PD_m14m98m4+HC_gr_PD_m52m58p24+HC_gr_PD_p22p2p10, data=MeanUptake)
MeanUptake$pred_probs_hypo<-predict(logit_model_hypo,type="response")
#increased clusters
logit_model_hyper<-glm(Fall ~ HC_kl_PD_p5m28p61+HC_kl_PD_p35m2p3+HC_kl_PD_m33m9p4, data=MeanUptake)
MeanUptake$pred_probs_hyper<-predict(logit_model_hyper,type="response")
#define object to plot
rocobj <- roc(MeanUptake$Fall,MeanUptake$pred_probs)## Setting levels: control = 0, case = 1
## Setting direction: controls < casesrocobj2<-roc(MeanUptake$Fall,MeanUptake$pred_probs_hypo)## Setting levels: control = 0, case = 1
## Setting direction: controls < casesrocobj3<-roc(MeanUptake$Fall,MeanUptake$pred_probs_hyper)## Setting levels: control = 0, case = 1
## Setting direction: controls < casesauc <- round(auc(MeanUptake$Fall, MeanUptake$pred_probs),4) ## Setting levels: control = 0, case = 1
## Setting direction: controls < casesauc1 <- round(auc(MeanUptake$Fall, MeanUptake$pred_probs_hypo),4) ## Setting levels: control = 0, case = 1
## Setting direction: controls < casesauc2 <- round(auc(MeanUptake$Fall, MeanUptake$pred_probs_hyper),4) ## Setting levels: control = 0, case = 1
## Setting direction: controls < casesci_roc1<-ci.auc(rocobj)
ci_roc2<-ci.auc(rocobj2)
ci_roc3<-ci.auc(rocobj3)
ci_band1<-ci.se(rocobj, specificities=seq(0,1,length.out=100))
ci_band2<-ci.se(rocobj2, specificities=seq(0,1,length.out=100))
ci_band3<-ci.se(rocobj3, specificities=seq(0,1,length.out=100))
ci_data<-data.frame(specificity=as.numeric(rownames(ci_band1)),
sensitivity=ci_band1[,2],
lower=ci_band1[,1],
upper=ci_band1[,3])
ci_data2<-data.frame(specificity=as.numeric(rownames(ci_band2)),
sensitivity=ci_band2[,2],
lower=ci_band2[,1],
upper=ci_band2[,3])
ci_data3<-data.frame(specificity=as.numeric(rownames(ci_band3)),
sensitivity=ci_band3[,2],
lower=ci_band3[,1],
upper=ci_band3[,3])
#create ROC plot
ggroc(list("All clusters" = rocobj,
"Hypometabolic" = rocobj2,
"Hypermetabolic" = rocobj3), size = 1) +
ggtitle(paste0('ROC Curve ', '(AUC = ', auc, ')'))+
theme_minimal()+
geom_ribbon(data=ci_data,aes(x=specificity, ymin=lower, ymax=upper),
inherit.aes= FALSE,fill="grey", alpha=0.2)+
geom_ribbon(data=ci_data2,aes(x=specificity, ymin=lower, ymax=upper),
inherit.aes= FALSE, fill= "#1F78B4", alpha=0.2)+
geom_ribbon(data=ci_data3,aes(x=specificity, ymin=lower, ymax=upper),
inherit.aes= FALSE,fill= "#BFD200", alpha=0.2)+
scale_colour_manual(values = c("All clusters" = "gray",
"Hypometabolic" = "#1F78B4",
"Hypermetabolic" = "#BFD200"),
guide = guide_legend(title = NULL))+
theme(text = element_text(size = 30),
axis.text.x = element_text(size=20),
axis.text.y = element_text(size=20))+
theme(legend.position = c(0.6, 0.4))
ggsave("roc_allClusters_ci.png")## Saving 7 x 5 in imagesetwd("~/Documents/RSfPETCog")
#load excel file
RSfPET<-read.xlsx(xlsxFile = "RS_fPET_fMRT_Data.xlsx",colNames = TRUE)
#define matrix for variation coefficient
Vk<-matrix(ncol=27, nrow=12)
#loop over subjects to fill matrix with variataion coefficient per subject
#calculate variation in subjects
Subjects<-c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029")
for(i in 1:length(Subjects)){
Vk[1,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m10m16m16)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m10m16m16)*100
Vk[2,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m14m98m4)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m14m98m4)*100
Vk[3,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m16_m2_p12)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m16_m2_p12)*100
Vk[4,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m34_p48_p2)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m34_p48_p2)*100
Vk[5,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m52m58p24)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_m52m58p24)*100
Vk[6,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p14m96m16)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p14m96m16)*100
Vk[7,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p22p2p10)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p22p2p10
)*100
Vk[8,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p46m62p36)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p46m62p36)*100
Vk[9,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p8m16m16)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_gr_PD_p8m16m16)*100
Vk[10,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_m33m9p4)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_m33m9p4)*100
Vk[11,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_p5m28p61)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_p5m28p61)*100
Vk[12,i]<- sd(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_p35m2p3)/mean(subset(gmFlexData,record_id==Subjects[i])$HC_kl_PD_p35m2p3)*100
}
#transform into dataframe
Vk_data<-data.frame(
Vk)
colnames(Vk_data)<-Subjects
tVk_data<-t(Vk_data)
colnames(tVk_data)<-colnames(gmFlexData[,1:12])
tVk_data## HC_gr_PD_m10m16m16 HC_gr_PD_m14m98m4 HC_gr_PD_m16_m2_p12
## RS_fPET001 9.149458 1.609397 5.520428
## RS_fPET002 11.174374 2.240504 8.320740
## RS_fPET003 8.022620 1.932647 9.629352
## RS_fPET004 10.962013 2.140445 7.782423
## RS_fPET005 8.504113 1.658767 7.300665
## RS_fPET006 8.744335 2.024273 8.351591
## RS_fPET007 8.111821 1.420442 4.606805
## RS_fPET008 11.383303 1.591634 8.671552
## RS_fPET010 9.868059 1.547594 7.268533
## RS_fPET011 5.203967 1.281799 5.507259
## RS_fPET012 10.621481 2.419876 7.434115
## RS_fPET014 9.893266 1.621636 7.851973
## RS_fPET015 7.897227 1.797752 6.324236
## RS_fPET016 11.171252 1.820831 10.433023
## RS_fPET017 6.924190 1.404396 6.623234
## RS_fPET018 8.918926 1.798407 7.407190
## RS_fPET019 8.683440 1.533015 6.715557
## RS_fPET020 6.211866 1.078035 4.450474
## RS_fPET021 8.110944 1.420442 4.606790
## RS_fPET022 8.671518 1.801752 5.561755
## RS_fPET023 6.272215 1.597641 5.153214
## RS_fPET024 14.587214 2.313922 9.904105
## RS_fPET025 8.964104 1.921321 7.927108
## RS_fPET026 9.633082 1.345150 7.113614
## RS_fPET027 7.512787 1.529448 8.411124
## RS_fPET028 7.130410 1.565766 8.438449
## RS_fPET029 7.721407 1.761468 7.009617
## HC_gr_PD_m34_p48_p2 HC_gr_PD_m52m58p24 HC_gr_PD_p14m96m16
## RS_fPET001 2.116563 1.526527 2.986196
## RS_fPET002 2.339266 2.032020 2.874476
## RS_fPET003 1.428722 1.185980 2.478086
## RS_fPET004 2.413514 2.026804 2.747820
## RS_fPET005 2.327577 1.748098 2.923575
## RS_fPET006 2.092383 1.605132 2.224615
## RS_fPET007 1.266032 1.504451 1.456026
## RS_fPET008 2.676008 2.079021 2.604002
## RS_fPET010 2.194953 1.429761 2.496667
## RS_fPET011 2.107382 1.585552 1.429364
## RS_fPET012 6.745391 1.807977 2.899212
## RS_fPET014 1.734242 2.030424 3.190650
## RS_fPET015 1.974489 1.594484 2.193591
## RS_fPET016 2.329207 1.556435 3.368566
## RS_fPET017 1.796029 1.329104 1.613693
## RS_fPET018 1.611327 1.446840 2.177380
## RS_fPET019 1.571010 1.314334 1.852377
## RS_fPET020 1.162451 1.195506 1.225862
## RS_fPET021 1.265957 1.504512 1.456278
## RS_fPET022 1.672039 1.506442 1.868991
## RS_fPET023 1.973596 1.342139 2.538681
## RS_fPET024 2.003060 2.303311 3.850872
## RS_fPET025 1.761571 1.397084 2.770462
## RS_fPET026 1.919996 1.389658 2.613000
## RS_fPET027 2.469319 2.107585 2.203953
## RS_fPET028 3.513236 1.443132 2.615536
## RS_fPET029 1.713690 1.702189 2.463870
## HC_gr_PD_p22p2p10 HC_gr_PD_p46m62p36 HC_gr_PD_p8m16m16
## RS_fPET001 5.171534 2.312951 7.675469
## RS_fPET002 6.612009 2.694771 9.149929
## RS_fPET003 5.026277 4.465347 9.592589
## RS_fPET004 6.786384 2.613275 9.908314
## RS_fPET005 5.199675 2.624397 12.276443
## RS_fPET006 6.292950 1.994300 9.857803
## RS_fPET007 4.210123 1.764034 5.636190
## RS_fPET008 5.625373 2.548775 8.576520
## RS_fPET010 5.609316 2.323284 11.026218
## RS_fPET011 3.923005 1.291514 6.123950
## RS_fPET012 6.284234 3.037558 8.272361
## RS_fPET014 6.907067 2.370053 10.131723
## RS_fPET015 5.269210 1.876998 9.049145
## RS_fPET016 6.026022 2.199640 9.985203
## RS_fPET017 4.178573 2.228051 6.765833
## RS_fPET018 4.344602 2.128189 9.290340
## RS_fPET019 4.015049 1.417233 8.160187
## RS_fPET020 4.014575 1.301726 7.246938
## RS_fPET021 4.210229 1.763910 5.636299
## RS_fPET022 4.971593 2.474279 8.798017
## RS_fPET023 4.479406 1.398938 10.423135
## RS_fPET024 8.101545 3.915617 15.915154
## RS_fPET025 4.355873 2.249852 11.326078
## RS_fPET026 5.901364 2.265688 7.772430
## RS_fPET027 5.116952 1.716365 9.864777
## RS_fPET028 4.337098 2.322549 9.663267
## RS_fPET029 4.831891 1.989895 8.642143
## HC_kl_PD_m33m9p4 HC_kl_PD_p5m28p61 HC_kl_PD_p35m2p3
## RS_fPET001 6.605761 2.2677417 4.263425
## RS_fPET002 5.813709 2.0108309 4.140702
## RS_fPET003 6.637324 1.9997148 3.977227
## RS_fPET004 6.203789 1.5369858 4.011151
## RS_fPET005 5.917214 2.0229921 3.627611
## RS_fPET006 4.965329 1.3606797 2.920158
## RS_fPET007 4.702880 2.0652449 2.863745
## RS_fPET008 6.629785 2.2372681 3.851112
## RS_fPET010 4.722849 1.3858007 2.959791
## RS_fPET011 4.101901 0.8377055 2.757841
## RS_fPET012 8.266548 2.2360301 4.723090
## RS_fPET014 5.971227 2.5375296 3.942221
## RS_fPET015 5.650172 1.8496787 2.730661
## RS_fPET016 5.346292 1.5727147 4.449450
## RS_fPET017 4.666100 1.0056832 3.500288
## RS_fPET018 6.263513 1.7842069 4.216995
## RS_fPET019 5.991526 1.0462713 2.813564
## RS_fPET020 3.398688 0.9437292 2.754232
## RS_fPET021 4.703037 2.0654299 2.863897
## RS_fPET022 6.557544 1.3977717 4.405924
## RS_fPET023 5.684074 0.7907276 3.756343
## RS_fPET024 6.812369 1.6182024 7.303282
## RS_fPET025 5.584171 2.2873223 3.406044
## RS_fPET026 6.014291 1.1979635 3.045293
## RS_fPET027 6.803682 1.2740988 4.228673
## RS_fPET028 5.534481 1.3830695 2.675859
## RS_fPET029 4.566263 1.3834854 4.164847#create group variable
#in groups HC vs PD
ID_HCPD<-MeanUptake$diagnose
tVk_group<-as.data.frame(tVk_data)
tVk_group<-data.frame(lapply(tVk_group,as.numeric))
tVk_group$ID_HCPD<-ID_HCPD
tVk_group$ID_HCPD<-factor(tVk_group$ID_HCPD)
summary(tVk_group$ID_HCPD)## HC PD
## 13 14variables<- colnames(tVk_group[,1:12])
grp<-levels(tVk_group$ID_HCPD)
#function
perm_test<- function(variable, group_var){
group_var<-as.factor(group_var)
grp_levels<-levels(group_var)
obs_diff<-tapply(variable,group_var,mean, na.rm=TRUE)
obs_diff_A<-obs_diff[grp_levels[2]]-obs_diff[grp_levels[1]]
#perm
perm_diffsA<-replicate(10000,{
permuted_groups<-sample(group_var)
perm_mean_diff_A<-tapply(variable,permuted_groups, mean, na.rm=TRUE)
perm_mean_diff_A[grp_levels[2]]-perm_mean_diff_A[grp_levels[1]]
})
print(length(perm_diffsA))
#p-values
p_value_A<-mean(abs(perm_diffsA)>=abs(obs_diff_A))
return(c(p_value_A))
}
#save results
results<-data.frame(Variable=character(), p_value=numeric(), stringsAsFactors = FALSE)
#loop over variables
for (var in variables){
p_values<-perm_test(tVk_group[[var]], tVk_group$ID_HCPD)
results<-rbind(results, data.frame(Variable=var, p_value_A=p_values[1]))
}## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000
## [1] 10000print(results)## Variable p_value_A
## 1 HC_gr_PD_m10m16m16 0.3454
## 2 HC_gr_PD_m14m98m4 0.9045
## 3 HC_gr_PD_m16_m2_p12 0.2495
## 4 HC_gr_PD_m34_p48_p2 0.8238
## 5 HC_gr_PD_m52m58p24 0.8775
## 6 HC_gr_PD_p14m96m16 0.0760
## 7 HC_gr_PD_p22p2p10 0.6586
## 8 HC_gr_PD_p46m62p36 0.1359
## 9 HC_gr_PD_p8m16m16 0.7014
## 10 HC_kl_PD_m33m9p4 0.7980
## 11 HC_kl_PD_p5m28p61 0.0235
## 12 HC_kl_PD_p35m2p3 0.9347#gghalf plots vc
library(gghalves)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# For multiple boxplots we need additional long format
long <- tVk_group %>%
pivot_longer(
cols = colnames(gmFlexData[, 1:12]), # Assumes gmFlexData has these columns
names_to = "region",
values_to = "value"
)
# Calculate max value for each region to guide ordering
region_summary <- long %>%
group_by(region) %>%
summarise(max_y = max(value, na.rm = TRUE)) %>%
arrange(max_y) # Arrange by ascending max value (regions with lower max_y first)
# Use arrange(desc(max_y)) for descending order if preferred
# Get the ordered region names
ordered_region_names <- region_summary$region
# Re-level the region factor in the long dataframe
# This new order will be respected by facet_wrap
long$region <- factor(long$region, levels = ordered_region_names)
# Create the plot
ggplot(long, aes(x = region, y = value, fill = ID_HCPD)) +
geom_half_violin(side = "r", trim = FALSE, alpha = 0.5, show.legend = FALSE) +
geom_half_boxplot(side = "l") +
facet_wrap(~region,scales="free", # You can add ncol or nrow here e.g. ncol=4 for 4 columns
labeller = as_labeller(c(
"HC_gr_PD_m10m16m16" = "left SN", "HC_gr_PD_m14m98m4" = "left OC",
"HC_gr_PD_m16_m2_p12" = "left AP","HC_gr_PD_m34_p48_p2" ="left MFC",
"HC_gr_PD_m52m58p24"="left AG","HC_gr_PD_p14m96m16" ="right OC",
"HC_gr_PD_p22p2p10"="right AP","HC_gr_PD_p46m62p36"="right AG",
"HC_gr_PD_p8m16m16"="right SN" ,"HC_kl_PD_m33m9p4"="left PP",
"HC_kl_PD_p5m28p61"="MC", "HC_kl_PD_p35m2p3"="right PP"
)))+
scale_fill_manual(labels=c("HC","PD"),values=c("#1F78B4", "#BFD200"))+
theme_minimal()+
labs(title = "", x = "Region",
y = "VC")+
theme_minimal() +
labs(title = "", x = "Region", y = "VC") +
theme(text = element_text(size = 25),
title = element_text(size = 20),
axis.text.x = element_blank(),
axis.text.y = element_text(size=20),
legend.title = element_blank())
# save plot
ggsave("Vk_datadriven_boxplots.png", width = 12, height = 10, dpi = 300)
#plot Vk vs clinical scores
tVk_data<-as.data.frame(tVk_data)
tVk_group$ID_HCPD<-factor(ID_HCPD)
#vs MoCA in whole group
plot(tVk_group$HC_gr_PD_p14m96m16,subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum)
cor.test(tVk_group$HC_gr_PD_p14m96m16,subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum, method="spearman")## Warning in cor.test.default(tVk_group$HC_gr_PD_p14m96m16, subset(RSfPET, :
## Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: tVk_group$HC_gr_PD_p14m96m16 and subset(RSfPET, record_id != "RS_fPET013" & record_id != "RS_fPET009")$MoCa_Sum
## S = 4615, p-value = 0.03428
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.408731# Spearman's rank correlation rho
#data: tVk_data$HC_gr_PD_p14m96m16 and subset(RSfPET, record_id != "RS_fPET013" & record_id != "RS_fPET009")$MoCa_Sum
#S = 4615, p-value = 0.03428
#alternative hypothesis: true rho is not equal to 0
#sample estimates:
# rho
#-0.408731
model<-lm(subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum ~ tVk_group$HC_gr_PD_p14m96m16)
abline(model)
summary(model)##
## Call:
## lm(formula = subset(RSfPET, record_id != "RS_fPET013" & record_id !=
## "RS_fPET009")$MoCa_Sum ~ tVk_group$HC_gr_PD_p14m96m16)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9036 -1.4685 0.6058 2.4538 4.2903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.038 2.616 12.627 2.39e-12 ***
## tVk_group$HC_gr_PD_p14m96m16 -2.550 1.050 -2.429 0.0227 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.426 on 25 degrees of freedom
## Multiple R-squared: 0.1909, Adjusted R-squared: 0.1585
## F-statistic: 5.899 on 1 and 25 DF, p-value: 0.02268#MC cluster
cor.test(tVk_group$HC_kl_PD_p5m28p61,subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum, method="spearman")## Warning in cor.test.default(tVk_group$HC_kl_PD_p5m28p61, subset(RSfPET, :
## Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: tVk_group$HC_kl_PD_p5m28p61 and subset(RSfPET, record_id != "RS_fPET013" & record_id != "RS_fPET009")$MoCa_Sum
## S = 4716.1, p-value = 0.02177
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4396019# Spearman's rank correlation rho
#data: tVk_data$HC_kl_PD_p5m28p61 and subset(RSfPET, record_id != "RS_fPET013" & record_id != "RS_fPET009")$MoCa_Sum
#S = 4716.1, p-value = 0.02177
#alternative hypothesis: true rho is not equal to 0
#sample estimates:
# rho
#-0.4396019
#with ggplot
ggplot(data=tVk_data,aes(x=HC_kl_PD_p5m28p61, y=subset(RSfPET,record_id!="RS_fPET013"&record_id!="RS_fPET009")$MoCa_Sum, color=ID_HCPD)) +
geom_point()+
geom_smooth(method='lm', formula= y~x,show.legend = FALSE)+
scale_color_manual(labels=c("HC","PD"),values=c("#1F78B4", "#BFD200"))+
theme_bw()+
theme(legend.position = c(0.9, 0.8), legend.title=element_blank())+
labs(title = "", x = "Motor Cortex VC",
y = "MoCa sum score")+
theme(text = element_text(size = 30),
axis.text.x = element_text(size=20),
axis.text.y = element_text(size=20))
#with Cog Z from demographic script
load("RSfPET.RData")
cor.test(subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z,tVk_group$HC_kl_PD_p5m28p61, method="spearman")##
## Spearman's rank correlation rho
##
## data: subset(RSfPET, record_id %in% c("RS_fPET001", "RS_fPET002", "RS_fPET003", "RS_fPET004", "RS_fPET005", "RS_fPET006", "RS_fPET007", "RS_fPET008", "RS_fPET010", "RS_fPET011", "RS_fPET012", "RS_fPET014", "RS_fPET015", "RS_fPET016", "RS_fPET017", "RS_fPET018", "RS_fPET019", "RS_fPET020", "RS_fPET021", "RS_fPET022", "RS_fPET023", "RS_fPET024", "RS_fPET025", "RS_fPET026", "RS_fPET027", "RS_fPET028", "RS_fPET029"))$Cog_z and tVk_group$HC_kl_PD_p5m28p61
## S = 4700, p-value = 0.02443
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4346764#t = -2.7353, df = 25, p-value = 0.01129
#alternative hypothesis: true correlation is not equal to 0
#95 percent confidence interval:
# -0.7273049 -0.1222158
#sample estimates:
# cor
#-0.4799396
#Figure
#with ggplot
ggplot(data=tVk_group,aes(x=HC_kl_PD_p5m28p61, y=subset(RSfPET,record_id%in%c("RS_fPET001","RS_fPET002","RS_fPET003","RS_fPET004","RS_fPET005","RS_fPET006","RS_fPET007","RS_fPET008","RS_fPET010","RS_fPET011","RS_fPET012","RS_fPET014","RS_fPET015","RS_fPET016","RS_fPET017","RS_fPET018","RS_fPET019","RS_fPET020","RS_fPET021","RS_fPET022","RS_fPET023","RS_fPET024","RS_fPET025","RS_fPET026","RS_fPET027","RS_fPET028","RS_fPET029"))$Cog_z,tVk_data$HC_kl_PD_p5m28p61, color=ID_HCPD, shape=RSfPET$MCI_NC_HC)) +
geom_point(size=3)+
geom_smooth(data=tVk_data,aes(group=1),method='lm', formula= y~x,show.legend = FALSE, color="black", se=FALSE)+
scale_color_manual(labels=c("HC","PD"),values=c("#1F78B4", "#BFD200"))+
theme_minimal()+
# theme(legend.position = "none", #legend.title=element_blank())+
labs(title = "", x = "Motor Cortex VC",
y = "Cognition z-score",
color="Group",
shape="Group")+
theme(text = element_text(size = 30),
axis.text.x = element_text(size=20),
axis.text.y = element_text(size=20))+
scale_shape_manual(values=c("HC"=16, "MCI"=17, "PD-NC"=16,"PD-MCI"=17))+
geom_text(x=1.5,y=-1.0, label="r = -0.48, p = 0.011", color="black", size=7)## Warning: The following aesthetics were dropped during statistical transformation: shape.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("Legend.png")## Saving 7 x 5 in image## Warning: The following aesthetics were dropped during statistical transformation: shape.
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).